Notebook demonstrating how to load data from the OpenScope predictive coding project¶
This project collected 2-photon calcium imaging data in 3 cortical areas (VISp, VISpm, RSP) across multiple depths in response to several stimulus types, including predictable image sequences with unexpected oddball images, spatial occlusion of varying degrees, natural movies, and control conditions with shuffled sequence and oddball images. In each mouse, a retrograde tracer was injected in either VISp or RSP to enable identification of feedforward and feedback projection neurons between areas. Biometric data including mouse running behavior and pupil diameter were also acquired.
In [1]:
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_context('notebook', font_scale=1.5, rc={'lines.markeredgewidth': 2})
sns.set_style('white', {'axes.spines.right': False, 'axes.spines.top': False, 'xtick.bottom': False, 'ytick.left': False,})
In [2]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
select an experiment from the manifest¶
In [3]:
manifest_file = r"C:\Users\marinag\Dropbox\opc_analysis\opc_production_manifest.xlsx"
manifest = pd.read_excel(manifest_file)
In [4]:
# limit to experiments that passed QC
manifest = manifest[manifest['experiment_state']=='passed']
In [5]:
print('total mice:', len(manifest['mouse_id'].unique()))
print('# V1 injections:', len(manifest[manifest['injection_area']=='VISp']['mouse_id'].unique()))
print('# RSP injections:', len(manifest[manifest['injection_area']=='RSP']['mouse_id'].unique()))
print('# V1 imaging sessions:', len(manifest[manifest['imaging_area']=='VISp']['experiment_id'].unique()))
print('# RSP imaging sessions:', len(manifest[manifest['imaging_area']=='RSP']['experiment_id'].unique()))
print('# PM imaging sessions:', len(manifest[manifest['imaging_area']=='VISpm']['experiment_id'].unique()))
In [6]:
VISpm_experiments = manifest[(manifest.imaging_area=='VISpm')].experiment_id.values
experiment_id = int(VISpm_experiments[5])
In [7]:
# info for this experiment
manifest[manifest.experiment_id==experiment_id]
Out[7]:
load data for a single 2-photon imaging session¶
In [8]:
# data streams for each imaging session are saved as .h5 files in each experiments's folder in this directory
cache_dir = r'C:\Users\marinag\Dropbox\opc_analysis'
In [9]:
# the OpenScopePredictiveCodingDataset class loads the relevant .h5 files and sets the data as attributes
from openscope_predictive_coding.ophys.dataset.openscope_predictive_coding_dataset import OpenScopePredictiveCodingDataset
dataset = OpenScopePredictiveCodingDataset(experiment_id, cache_dir)
In [10]:
from openscope_predictive_coding.ophys.response_analysis.response_analysis import ResponseAnalysis
# set use_events to True if you want responses in events instead of dF/F
analysis = ResponseAnalysis(dataset, use_events=True)
In [11]:
# plot max intensity projection for this experiment
plt.imshow(dataset.max_projection, cmap='gray')
Out[11]:
In [12]:
# get cell traces
dataset.dff_traces.head(5)
Out[12]:
In [13]:
# plot a heatmap of all cell traces (this is slow)
fig, ax = plt.subplots(figsize=(20,5))
sns.heatmap(dataset.dff_traces_array, vmax=np.percentile(dataset.dff_traces_array, 99), ax=ax, linecolor=None, linewidth=0)
ax.set_ylabel('cell_index')
ax.set_xlabel('2P frames')
Out[13]:
In [14]:
def plot_sorted_traces_heatmap(dataset, analysis, ax=None, save=False, use_events=False):
import openscope_predictive_coding.ophys.response_analysis.response_processing as rp
import openscope_predictive_coding.ophys.plotting.experiment_summary_figures as esf
if use_events:
traces = dataset.events_array.copy()
traces = rp.filter_events_array(traces, scale=2)
traces = esf.reorder_traces(traces, analysis)
vmax = 0.01
# vmax = np.percentile(traces, 99)
label = 'event magnitude'
suffix = '_events'
else:
traces = dataset.dff_traces_array
traces = esf.reorder_traces(traces, analysis)
vmax = np.percentile(traces, 99)
label = 'dF/F'
suffix = ''
if ax is None:
figsize = (14, 5)
fig, ax = plt.subplots(figsize=figsize)
cax = ax.pcolormesh(traces, cmap='magma', vmin=0, vmax=vmax)
ax.set_ylabel('cells')
interval_seconds = 5 * 60
ophys_frame_rate = int(dataset.metadata.ophys_frame_rate.values[0])
upper_limit, time_interval, frame_interval = esf.get_upper_limit_and_intervals(traces, dataset.ophys_timestamps,
ophys_frame_rate)
ax.set_xticks(np.arange(0, upper_limit, interval_seconds * ophys_frame_rate))
ax.set_xticklabels(np.arange(0, upper_limit / ophys_frame_rate, interval_seconds))
ax.set_xlabel('time (seconds)')
cb = plt.colorbar(cax, pad=0.015)
cb.set_label(label, labelpad=3)
if save:
esf.save_figure(fig, figsize, dataset.analysis_dir, 'experiment_summary',
str(dataset.experiment_id) + 'sorted_traces_heatmap' + suffix)
return ax
In [ ]:
In [15]:
# plot_sorted_traces_heatmap(dataset, analysis, ax=None, save=False, use_events=False)
rastermap¶
In [16]:
fig, ax = plt.subplots(figsize=(20,5))
sns.heatmap(dataset.dff_traces_array, vmax=np.percentile(dataset.dff_traces_array, 99), ax=ax, linecolor=None, linewidth=0)
Out[16]:
In [17]:
from rastermap import Rastermap
model = Rastermap(n_components=2, n_X=10, nPC=100, init='pca')
# fit does not return anything, it adds attributes to model
# attributes: embedding, u, s, v, isort1
sp = dataset.dff_traces_array.copy()
model.fit(sp)
plt.imshow(sp[model.isort, :])
# fit_transform returns embedding (upsampled cluster identities)
embedding = model.fit_transform(sp)
# transform can be used on new samples with the same number of features as sp
# embed2 = model.transform(sp2)
In [18]:
fig, ax = plt.subplots(figsize=(20,5))
sns.heatmap(sp[model.isort, :], vmax=np.percentile(dataset.dff_traces_array, 99), ax=ax, linecolor=None, linewidth=0)
Out[18]:
In [19]:
sns.heatmap(embedding)
Out[19]:
across experiments¶
In [20]:
cache_dir = r'\\allen\programs\braintv\workgroups\nc-ophys\opc\opc_analysis'
In [21]:
from rastermap import Rastermap
model = Rastermap(n_components=2, n_X=10, nPC=100, init='pca')
# fit does not return anything, it adds attributes to model
# attributes: embedding, u, s, v, isort1
for experiment_id in manifest.experiment_id.unique():
experiment_id = int(experiment_id)
try:
dataset = OpenScopePredictiveCodingDataset(experiment_id, cache_dir)
sp = dataset.dff_traces_array.copy()
model.fit(sp)
# plt.imshow(sp[model.isort, :])
# fit_transform returns embedding (upsampled cluster identities)
embedding = model.fit_transform(sp)
# transform can be used on new samples with the same number of features as sp
# embed2 = model.transform(sp2)
fig, ax = plt.subplots(figsize=(20,5))
ax = sns.heatmap(sp[model.isort, :], vmax=np.percentile(dataset.dff_traces_array, 99), ax=ax, linecolor=None, linewidth=0)
ax.set_title(dataset.analysis_folder)
except:
print('problem for',experiment_id)
In [22]:
model.isort
Out[22]:
In [ ]:
dataset.metadata
In [ ]:
In [ ]:
In [ ]:
In [ ]:
# minmax standardized
In [ ]:
methods = ['single', 'average', 'weighted', 'median', 'ward']
metrics = ['euclidean', 'correlation', 'mahalanobis']
In [ ]:
sns.clustermap(dataset.dff_traces_array, method='average', metric='correlation', z_score=None, standard_scale=None)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
for method in methods:
for metric in metrics:
try:
fig, ax = plt.subplots()
ax = sns.clustermap(dataset.dff_traces_array, method='average', metric='euclidean', z_score=None, standard_scale=0, ax=ax)
ax.set_title(method+', '+metric+', standard scale')
except:
pass
In [ ]:
sns.clustermap(dataset.dff_traces_array, method='average', metric='euclidean', z_score=None, standard_scale=0,)
In [ ]:
sns.clustermap(dataset.dff_traces_array, method='average', metric='euclidean', z_score=0, standard_scale=None)
In [ ]:
In [ ]:
In [14]:
# plot one cell's dF/F trace with x-axis in time in seconds
cell_specimen_id = dataset.cell_specimen_ids[0]
dff_trace = dataset.dff_traces.loc[cell_specimen_id].values[0]
timestamps = dataset.timestamps_ophys
fig, ax = plt.subplots(figsize=(20,3))
ax.plot(timestamps, dff_trace)
ax.set_ylabel('dF/F')
ax.set_xlabel('time (sec)')
Out[14]:
In [15]:
# plot deconvolved calcium events for the same cell
events_trace = dataset.events.loc[cell_specimen_id].values[0]
timestamps = dataset.timestamps_ophys
fig, ax = plt.subplots(figsize=(20,3))
ax.plot(timestamps, events_trace)
ax.set_ylabel('dF/F')
ax.set_xlabel('time (sec)')
Out[15]:
In [16]:
# plot mouse running speed with x-axis in seconds
run_speed = dataset.running_speed.speed.values
timestamps = dataset.running_speed.time.values
fig, ax = plt.subplots(figsize=(20,3))
ax.plot(timestamps, run_speed)
ax.set_ylabel('running speed (cm/s)')
ax.set_xlabel('time (sec)')
Out[16]:
looks like this mouse didnt run
In [17]:
# what are the different stimulus types?
dataset.stimulus_table.session_block_name.unique()
Out[17]:
In [18]:
# plot a cell trace with duration of stimulus blocks indicated
dff_trace = dataset.dff_traces.loc[cell_specimen_id].values[0]
timestamps = dataset.timestamps_ophys
fig, ax = plt.subplots(figsize=(20,3))
ax.plot(timestamps, dff_trace)
ax.set_ylabel('dF/F')
ax.set_xlabel('time (sec)')
block_df = dataset.stimulus_block_table.copy()
colors = sns.color_palette('deep')
for i, block_name in enumerate(block_df.block_name.values):
start_time = block_df[block_df.block_name==block_name].start_time.values[0]
end_time = block_df[block_df.block_name==block_name].end_time.values[0]
ax.axvspan(start_time, end_time, facecolor=colors[i], edgecolor='none', alpha=0.3, linewidth=0, zorder=1, label=block_name)
ax.legend(bbox_to_anchor=(1,1))
Out[18]:
In [19]:
# stimulus table has times of all stimulus presentations for each of the stimulus types
dataset.stimulus_table.head(3)
Out[19]:
get cell responses for each stimulus presentation¶
the ResponseAnalysis class does the work of temporal alignment between 2P times and stimulus times and creates dataframes with the response of every cell for every stimulus in a given stimulus block type. If these dataframes have already been generated, the get_response_df() method loads the saved dataframe.
In [20]:
from openscope_predictive_coding.ophys.response_analysis.response_analysis import ResponseAnalysis
In [21]:
# set use_events to True if you want responses in events instead of dF/F
analysis = ResponseAnalysis(dataset, use_events=False)
In [22]:
# what are the stimulus types again?
dataset.stimulus_block_table
Out[22]:
In [23]:
# provide the block name to this method to get the stimulus triggered responses for that stimulus type
odf = analysis.get_response_df('oddball')
In [24]:
odf.head(3)
Out[24]:
In [25]:
odf.keys()
Out[25]:
In [26]:
# plot single trial response for one cell
trace = odf[(odf.cell_specimen_id==cell_specimen_id)].trace.values[50]
plt.plot(trace)
plt.xlabel('2P frames')
plt.ylabel('dF/F')
Out[26]:
In [27]:
# plot trial averaged response for all oddball images for one cell
mean_trace = odf[(odf.cell_specimen_id==cell_specimen_id)&(odf.oddball==True)].trace.mean()
times = odf[(odf.cell_specimen_id==cell_specimen_id)&(odf.oddball==True)].trace_timestamps.values[0]
fig, ax = plt.subplots()
ax.plot(times, mean_trace)
ax.set_xlabel('time after stimulus onset (sec)')
ax.set_ylabel('dF/F')
Out[27]:
In [28]:
# plot the populatin response vector for one trial - mean response of all cells to that stimulus
# stimulus_presentations_id is the index of each individual stimulus presentation across the entire session
oddball_trials = odf[odf.oddball==True].stimulus_presentations_id.unique()
mean_responses = odf[(odf.stimulus_presentations_id==oddball_trials[9])].mean_response.values
plt.plot(mean_responses, '.')
plt.xlabel('cells')
plt.ylabel('mean response')
plt.title('population response for one oddball trial')
Out[28]:
In [29]:
# how many oddball images are there anyway?
print(len(analysis.oddball_images), 'oddball images')
print('oddball image IDs are', analysis.oddball_images)
In [30]:
# what are the sequence images?
print('sequence image IDs are', analysis.sequence_images)
In [31]:
# plot one of the images
stimulus_template = dataset.get_stimulus_template()
image_id = analysis.sequence_images[1]
plt.imshow(stimulus_template[image_id])
Out[31]:
get trial averaged responses for a set of conditions¶
In [32]:
# start with dataframe with every cell's response to every stimulus presentation
odf = analysis.get_response_df('oddball')
In [33]:
# utilities has some useful functions for trial averaging, determining significance of response, etc.
import openscope_predictive_coding.ophys.response_analysis.utilities as ut
In [34]:
# the get_mean_df() function gets trial averaged responses and other metrics for a set of conditions
conditions=['cell_specimen_id', 'oddball', 'image_id'] # columns in response dataframe over which to groupby
oddball_mean_df = ut.get_mean_df(odf, conditions=conditions)
In [35]:
oddball_mean_df.head()
Out[35]:
In [36]:
# plot the population average response for oddball and sequence images
last_in_sequence = analysis.sequence_images[-1]
plt.plot(oddball_mean_df[oddball_mean_df.oddball==True].mean_trace.mean(), label='oddball')
plt.plot(oddball_mean_df[oddball_mean_df.image_id==last_in_sequence].mean_trace.mean(), label='sequence')
plt.title('population average response')
plt.xlabel('2P frames')
plt.ylabel('dF/F')
plt.legend()
Out[36]:
In [37]:
# plot mean response across oddball images for some cell
cell_data = oddball_mean_df[oddball_mean_df.cell_specimen_id==cell_specimen_id]
for image_id in analysis.oddball_images:
trace = cell_data[cell_data.image_id==image_id].mean_trace.values[0]
plt.plot(trace, label=int(image_id))
plt.xlabel('2P frames')
plt.ylabel('dF/F')
plt.title('cell '+str(cell_specimen_id)+'\nmean oddball image responses')
plt.legend(bbox_to_anchor=(1,1), fontsize='medium', title='image_id')
Out[37]:
In [ ]:
In [ ]: